#### prepare packages and functions ####
require(estimatr)
require(dichromat)

issue.var.name <- c("a.article9", "b.defense", "c.revisionism", "d.women", 
                    "e.gay", "f.immigrant", "g.growth", "h.tax")

blue.orange <- colorschemes$BluetoOrange.10
group.col.1 <- c(blue.orange[10], "gray30", blue.orange[1])
group.col.2 <- c(paste0(blue.orange[7], "80"), "#CCCCCC80", 
                 paste0(blue.orange[4], "80"))
group.lty <- c(1, 2, 4)

issue.label <- c("Revision of Article 9", "Increasing Defense Power", 
                 "Historical Revisionism", "Women's Empowerment", 
                 "Gay Marriage", "Accepting Foreign Workers", 
                 "Economic Growth", "Progressive Tax")

# function for the quadratic regression
MM.predict <- function(data, attribute, outcome, predictor, values, id) {
  lm.data <- data.frame(id = data[, id], A = data[, attribute],
                        Y = data[, outcome], X = data[, predictor])
  l <- nlevels(lm.data$A)
  N <- max(lm.data$id)
  lm.result <- array(NA, c(3, 3, l))
  prediction <- array(NA, c(length(values), l, 3))
  for (i in 1:l) {
    level.label <- levels(lm.data$A)[i]
    QR.result <- lm_robust(Y ~ X + I(X ^ 2), data = lm.data,
                           subset = A == level.label,
                           clusters = id)
    lm.result[, 1, i] <- QR.result$coefficients
    lm.result[, 2, i] <- QR.result$std.error
    lm.result[, 3, i] <- QR.result$p.value
    prediction[, i, ] <- predict(QR.result, data.frame(X = values),
                                 interval = "confidence")$fit
  }
  list(lm.result = lm.result, prediction = prediction)
}

# function to conduct F-tests for linear models
age.1.F.test <- function(variable, data, choice) {
  if (choice == TRUE) {
    lm.formula <- as.formula(paste0("choice ~ age + ", variable))
  } else {
    lm.formula <- as.formula(paste0("rating ~ age + ", variable))
  }
  lm.model.1 <- lm(lm.formula, data = data)
  if (choice == TRUE) {
    lm.formula <- as.formula(paste0("choice ~ age * ", variable))
  } else {
    lm.formula <- as.formula(paste0("rating ~ age * ", variable))
  }
  lm.model.2 <- lm(lm.formula, data = data)
  anova.result <- anova(lm.model.1, lm.model.2)
  c(anova.result$"F"[2], anova.result$"Pr(>F)"[2])
}

# function to conduct F-tests for quadratic models
age.2.F.test <- function(variable, data, choice) {
  if (choice == TRUE) {
    lm.formula <- as.formula(paste0("choice ~ age + I(age ^ 2) + ", variable))
  } else {
    lm.formula <- as.formula(paste0("rating ~ age + I(age ^ 2) + ", variable))
  }
  lm.model.1 <- lm(lm.formula, data = data)
  if (choice == TRUE) {
    lm.formula <- as.formula(paste0("choice ~ (age + I(age ^ 2)) * ", variable))
  } else {
    lm.formula <- as.formula(paste0("rating ~ (age + I(age ^ 2)) * ", variable))
  }
  lm.model.2 <- lm(lm.formula, data = data)
  anova.result <- anova(lm.model.1, lm.model.2)
  c(anova.result$"F"[2], anova.result$"Pr(>F)"[2])
}

#### load data ####
task.data <- read.csv("task_data.csv")

# reorder the levels of attribute variables
position.label <- c("Agree", "Neither", "Disagree")
task.data$a.article9 <- factor(task.data$a.article9, levels = position.label)
task.data$b.defense <- factor(task.data$b.defense, levels = position.label)
task.data$c.revisionism <- factor(task.data$c.revisionism, levels = position.label)
task.data$d.women <- factor(task.data$d.women, levels = position.label)
task.data$e.gay <- factor(task.data$e.gay, levels = position.label)
task.data$f.immigrant <- factor(task.data$f.immigrant, levels = position.label)
task.data$g.growth <- factor(task.data$g.growth, levels = position.label)
task.data$h.tax <- factor(task.data$h.tax, levels = position.label)

#### analysis ####
## predict MMs by age (divided by 10)
choice.predict <- rating.predict <- list()
for (i in 1:8) {
  choice.predict[[i]] <- MM.predict(task.data, issue.var.name[i], "choice",
                                    "age.decimal", c(18:69) / 10, 
                                    "respondent.id")
  rating.predict[[i]] <- MM.predict(task.data, issue.var.name[i], "rating",
                                    "age.decimal", c(18:69) / 10, 
                                    "respondent.id")
}

# Figure 2
png("Figure_2.png", width = 6.2, height = 3.3, units = "in", pointsize = 8, 
    res = 2000, type = "cairo")
layout(matrix(1:8, 2, 4, byrow = TRUE))
par(mar = c(3, 3, 3, 1), lwd = 0.5)
for (i in 1:8) {
  plot(NULL, NULL, type = "n", bty = "n", xlim = c(18, 69), ylim = c(4.7, 6.3), 
       main = issue.label[i], xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(h = 5.5, lty = 3, col = "gray30")
  for (j in 3:1) {
    polygon(c(18:69, 69:18), 
            c(rating.predict[[i]]$prediction[, j, 2], 
              rev(rating.predict[[i]]$prediction[, j, 3])), 
            border = NA, col = group.col.2[j])
  }
  for (j in 3:1) {
    lines(18:69, rating.predict[[i]]$prediction[, j, 1], 
          lty = group.lty[j], col = group.col.1[j])
  }
  if (i == 1) {
    text(69, 6.27, "Agree", pos = 2)
    text(69, 5.72, "Neither", pos = 2)
    text(69, 5.15, "Disagree", pos = 2)
  }
  axis(1, at = c(seq(18, 58, 10), 69), lwd = 0.5)
  axis(2, at = seq(4.7, 6.3, 0.4), lwd = 0.5)
  mtext("Age", side = 1, line = 2, cex = 0.7)
  mtext("Marginal Mean", side = 2, line = 2, cex = 0.7)
}
dev.off()

# Figure A.4
cairo_pdf("Figure_A4.pdf", width = 6.2, height = 3.3, pointsize = 8)
layout(matrix(1:8, 2, 4, byrow = TRUE))
par(mar = c(3, 3, 3, 1), lwd = 0.5)
for (i in 1:8) {
  plot(NULL, NULL, type = "n", bty = "n", xlim = c(18, 69), ylim = c(0.3, 0.7), 
       main = issue.label[i], xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  for (j in 3:1) {
    polygon(c(18:69, 69:18), 
            c(choice.predict[[i]]$prediction[, j, 2], 
              rev(choice.predict[[i]]$prediction[, j, 3])), 
            border = NA, col = group.col.2[j])
  }
  for (j in 3:1) {
    lines(18:69, choice.predict[[i]]$prediction[, j, 1], 
          lty = group.lty[j], col = group.col.1[j])
  }
  if (i == 1) {
    text(69, 0.61, "Agree", pos = 2)
    text(69, 0.51, "Neither", pos = 2)
    text(69, 0.38, "Disagree", pos = 2)
  }
  axis(1, at = c(seq(18, 58, 10), 69), lwd = 0.5)
  axis(2, at = seq(0.3, 0.7, 0.1), lwd = 0.5)
  mtext("Age", side = 1, line = 2, cex = 0.7)
  mtext("Marginal Mean", side = 2, line = 2, cex = 0.7)
}
dev.off()

# Table A.4
for (i in 1:8) {
  print(issue.var.name[i])
  print(round(cbind(choice.predict[[i]]$lm.result[, , 1], 
                    choice.predict[[i]]$lm.result[, , 2], 
                    choice.predict[[i]]$lm.result[, , 3]), 3))
}

# Table A.5
for (i in 1:8) {
  print(issue.var.name[i])
  print(round(cbind(rating.predict[[i]]$lm.result[, , 1], 
                    rating.predict[[i]]$lm.result[, , 2], 
                    rating.predict[[i]]$lm.result[, , 3]), 3))
}

## conduct F-tests
# linear models
age.1.F.test.p.choice <- age.1.F.test.p.rating <- matrix(NA, 8, 2)
for (i in 1:8) {
  age.1.F.test.p.choice[i, ] <- age.1.F.test(issue.var.name[i], 
                                             data = task.data, choice = TRUE)
  age.1.F.test.p.rating[i, ] <- age.1.F.test(issue.var.name[i], 
                                             data = task.data, choice = FALSE)
}
rownames(age.1.F.test.p.choice) <- 
  rownames(age.1.F.test.p.rating) <- issue.var.name

# Table A.6-A
round(age.1.F.test.p.choice, 3)
round(age.1.F.test.p.rating, 3)

# quadratic models
age.2.F.test.p.choice <- age.2.F.test.p.rating <- matrix(NA, 8, 2)
for (i in 1:8) {
  age.2.F.test.p.choice[i, ] <- age.2.F.test(issue.var.name[i], 
                                             data = task.data, choice = TRUE)
  age.2.F.test.p.rating[i, ] <- age.2.F.test(issue.var.name[i], 
                                             data = task.data, choice = FALSE)
}
rownames(age.2.F.test.p.choice) <- 
  rownames(age.2.F.test.p.rating) <- issue.var.name

# Table A.6-B
round(age.2.F.test.p.choice, 3)
round(age.2.F.test.p.rating, 3)